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ABSTRACT 

We present a comprehensive study for the gamma-ray luminosity function (GLF) of blazars and their con- 
tribution to the extragalactic diffuse gamma-ray background (EGRB). Radio and gamma-ray luminosity cor- 
relation is introduced with a modest dispersion consistent with observations, to take into account the radio 
detectability which is important for the blazar identification. Previous studies considered only pure luminosity 
evolution (PLE) or pure density evolution, but here we introduce the luminosity dependent density evolution 
(LDDE) model, which is favored from the evolution of X-ray luminosity function (XLF) of AGNs. The model 
parameters are constrained by likelihood analyses about the observed redshift and gamma-ray flux distributions 
of the EGRET blazars. Interestingly, we find that the LDDE model gives a better fit to the observed distribu- 
tions than the PLE model, indicating that the LDDE model is also appropriate for gamma-ray blazars, and 
that the jet activity is universally correlated with the accretion history of AGNs. The normalization between 
the GLF and XLF is consistent with the unified picture of AGNs, when the beaming and a reasonable duty 
cycle of jet activity are taken into account. We then find that only 25-50% of the EGRB can be explained by 
unresolved blazars with the best-fit LDDE parameters. Unresolved blazars can account for all the EGRB only 
with a steeper index of the faint-end slope of the GLF, which is marginally consistent with the EGRET data but 
inconsistent with that of the XLF. Therefore unresolved AGNs cannot be the dominant source of the EGRB, 
unless there is a new population of gamma-ray emitting AGNs that evolves differently from the XLF of AGNs. 
Predictions for the GLAST mission are made, and we find that the best-fit LDDE model predicts about 3000 
blazars in the entire sky, which is considerably fewer (by a factor of more than three) than a previous estimate. 
Subject headings: diffuse radiation — galaxies: active — galaxies: evolution — galaxies: luminosity function 
— gamma rays: theory — quasars: general 



1. INTRODUCTION 

The origin of the extragalactic diffuse gamma-ray back- 
ground (EGRB) is one of the unsolved problems in astro- 
physics. The EGRB was first discovered by the SAS 2 satellite 
(Fichtel, Simpson, & Thompson 1978; Thompson & Fichtel 
1982) and subsequently confirmed by the Energetic Gamma 
Ray Experiment Telescope (EGRET) instrument aboard the 
Compton Gamma Ray Observatory ( CGRO). In the first anal- 
ysis of the EGRET data, the flux of the EGRB integrated 
above 100 MeV was determined to be (1.45 ±0.05) x 10"^ 
photons cm^^ s^' sr^' (Sreekumar et al. 1998). However, 
this value is strongly dependent on the modeling of the Galac- 
tic background, which is dominated by cosmic -ray interac- 
tions in the Galactic disk and must be subtracted from the 
background data. The latest analysis, which used a new model 
of the Galactic background, resulted in a slightly smaller 
value of the EGRB, (1.14±0.12) x 10"^ photons cm 
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sr"^ (Strong et al. 2004). 

EGRET detected many extragalactic high-energy gamma- 
ray sources that have been identified as active galactic nuclei 
(AGNs). Most of them fall into the blazar class of AGNs, 
and this is the only one extragalactic population confirmed in 
the third EGRET catalog (Hartman et al. 1999), constituting 
about 15% of the EGRB flux. Therefore unresolved blazars 
are the most likely candidate for the origin of the EGRB, at 
least contributing substantially, and this issue has been stud- 
ied in a number of papers (Padovani et al. 1993; Stecker, 
Salamon, & Malkan 1993; Salamon & Stecker 1994; Chiang 
et al. 1995; Stecker & Salamon 1996; Chiang & Mukherjee 



1998; Miicke & Pohl 2000). On the other hand, several alter- 
native candidates for the EGRB components have been pro- 
posed, e.g., intergalactic shocks produced by the formation of 
large-scale cosmological structures (Loeb & Waxman 2000; 
Totani & Kitayama 2000; Miniati 2002; Scharf & Mukherjee 
2002; Gabici & Blasi 2003; Keshet et al. 2003), or dark matter 
annihilation (Oda, Totani, & Nagashima 2005 and references 
therein). Therefore, it is important to determine whether the 
number of unresolved blazars are enough to account for all of 
the EGRB, but conclusions derived by the earlier studies are 
somewhat controversial. 

Stecker & Salamon (1996, hereafter SS96) estimated the 
unresolved blazar contribution with basic assumptions that 
EGRET blazars are the same population with flat-spectrum 
radio-loud quasars (FSRQs), and that the gamma-ray and ra- 
dio luminosities are linearly related. Then they constructed 
the blazar gamma-ray luminosity function (GLF) model from 
the FSRQ radio luminosity function (RLE), and found that 
blazars can account for 100% of the EGRB. However, their 
model was not compared with the available redshift distri- 
bution of the EGRET blazars, and hence it was uncertain 
whether this GLF model is statistically consistent with the 
EGRET blazar data. 

Quantitative comparison of GLF models to the flux and red- 
shift distributions of the EGRET blazars was performed by 
Chiang & Mukherjee (1998, hereafter CM98), and indeed, 
they found that the model of SS96 seriously overpredicts the 
number of low-redshift blazars detectable by the EGRET. 
CM98 then concluded that blazars can account for only 25% 
of the EGRB, based on the GLF model consistent with the 
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EGRET blazar distributions. ' 

However, the analysis of the GLF is not straightforward; 
a source of uncertainty is the detectability in the radio band. 
Most of the EGRET blazars are identified by finding radio 
counterparts, and hence they would remain unidentified if 
their radio counterparts are under the flux limit of radio sur- 
veys, even though their gamma-ray flux is above the EGRET 
sensitivity limit. Therefore one must estimate the probabil- 
ity of a model blazar having flux greater than the sensitivity 
limits, not only in the gamma-ray band but also in the radio 
band. CM98 introduced this probability in their analysis, but 
they assumed that there is no correlation between gamma-ray 
and radio luminosities of blazars. However, the assumption of 
no correlation at all over a wide range of gamma-ray and ra- 
dio luminosities induces some inconsistencies (see discussion 
given in Stecker & Salamon 2001), and it is physically reason- 
able to expect some level of correlation from the viewpoint of 
the standard synchrotron-inverse Compton model of blazars. 
Therefore we adopt a new treatment on this issue introducing 
a reasonable correlation (see ^2.3l for details). 

Miicke & Pohl (2000) approached this issue from the view- 
point of the unification scheme of radio-loud AGNs, which 
proposed that blazars are the beamed subclass of Fanaroff- 
Riley (PR) radio galaxies. They considered the GLF models 
based on the RLE of FR galaxies, and the flux and redshift 
distributions of blazars were used to constrain their GLF mod- 
els. Then they concluded that unresolved blazar contributuion 
to the EGRB is 20-40% assuming that the blazars extend to 
the maximum cutoff redshift of Zmax = 3, and 40-80% for 
Zmax — 5. However, the identification probability of a blazar, 
which may affect the estimate of the blazar contribution to the 
EGRB, was not incorporated in their analysis. 

To resolve this rather controversial situation, in this paper 
we make a comprehensive study of GLF models that are sta- 
tistically compared with the observed redshift and flux distri- 
butions, taking into account a reasonable correlation between 
gamma-ray and radio flux in a consistent manner with the ob- 
served gamma-ray to radio flux ratios. Then we make an es- 
timate of the blazar contribution to the EGRB flux, and also 
estimate the expected number of "unidentified blazars" due 
to the lack of radio detection, which can be compared with 
the number of high-Galactic-latitude unidentified sources in 
the third EGRET catalog. We also make some predictions for 
the future Gamma Ray Large Area Space Telescope (GLAST) 
observation, and discuss its prospects. 

In addition to these new aspects, we also try a new type 
of the GLF evolution model. The earlier studies treated the 
cosmological evolution of the blazar GLF as a pure luminos- 
ity evolution (PLE) or a pure density evolution. On the other 
hand, the cosmological evolution of the luminosity function 
of AGNs has been investigated intensively in soft X-ray (e.g., 
Miyaji, Hasinger, & Schmidt 2000, 2001; Hasinger, Miyaji, 
& Schmidt 2005) and hard X-ray (e.g., Boyle et al. 1998; La 
Franca et al. 2002; Cowie et al. 2003; Ueda et al. 2003) 
bands, recently. The former studies are mostly for type 1 
AGNs, while the latter ones are for both type 1 and 2 AGNs, 
since soft X-ray emission of type 2 AGNs is significantly ab- 
sorbed. These studies revealed that the overall behavior of 
the soft X-ray luminosity function (SXLF) and hard X-ray 

' In this paper, we refer the fraction of the blazar contribution in the EGRB 
always removing the already detected EGRET blazars, i.e., the ratio of the 
background flux from blazars under the EGRET detection limit to the EGRB 
flux of 1.14 X 10"^ photons cm"^ s"' reported by Strong et al. (2004). 



luminosity function (HXLF) of AGNs are very similar and 
best described with a luminosity dependent density evolu- 
tion (LDDE) where the peak redshift of density evolution in- 
creases with AGN luminosity (Miyaji, Hasinger, & Schmidt 
2000, 2001; Ueda et al. 2003; Hasinger, Miyaji, & Schmidt 
2005). Therefore it is reasonable to expect that the cosmolog- 
ical evolution of the blazar GLF may also be expressed by the 
LDDE. In this paper, we try two kinds of blazar GLF model; 
one is based on the FSRQ RLE (PLE model) and the other is 
based on the AGN SXLF (LDDE model). 

The paper will be organized as follows: In §2, we describe 
sample definition and formulations for the statistical analy- 
ses. In §3, we present models of the blazar GLF and results 
of comparison to the EGRET data. In §4, we address the 
prospects for the GLAST mission. Discussions and conclu- 
sions are given in §5 and §6, respectively. Throughout this 
paper, we adopt a ACDM universe with the density param- 
eter 51o = 0.3, the cosmological constant iljy ~ 0.7, and the 
Hubble constant Hq = 70 km s"' Mpc"-\ 

2. FORMULATIONS 
2.1. Sample Definition 

Most of blazars show considerable variability, but it is diffi- 
cult to incorporate in the statistical analysis of the GLF. There- 
fore we use the mean flux shown in the third EGRET cata- 
log. There are some blazars in the EGRET catalog whose sig- 
nificance of detection is above the EGRET threshold {4a for 
1^1 > 10°, and 5ct for \b\ < 10°) only in some viewing periods, 
and their mean significance of detection is under the EGRET 
threshold. We exclude such blazars from our analysis, and use 
only blazars whose mean significance of detection (as shown 
in the catalog) is above the EGRET threshold. There are 46 
blazars meeting this selection. 

Since the detection limit is sensitively dependent on the lo- 
cation in the sky, we will take it into account in the analysis. In 
the left panel of Figure ^ we plot the EGRET detection limit 
for each EGRET point source, which is calculated from mean 
flux and mean statistical significance of the detection by using 
equation (1) of CM98, against its Galactic latitude {crosses). 
The best-fit relation is also shown {solid line). Because of the 
variability, this fit may be different from the effective detec- 
tion threshold for blazars that we have chosen. In order to 
check this point, in the right panel of Figure ^ we plot the 
mean flux of the EGRET blazars, comparing them to the ob- 
tained detection limit curve. We see from this figure that most 
of blazars are above the detection limit curve, and hence it is 
reasonable to take this curve as the EGRET detection limit for 
the 46 blazars used in this paper. 

It is difficult to take into account the variety of blazar spec- 
tra, and we assume a single universal power-law spectrum 
for all blazars, with the photon spectral index of a = 2.2 
{dN^/dE.y o<: E^"). This value is close to the mean of the 
EGRETblazars(a = 2.15±0.04,Sreekumai-etal. 1998). We 
confirmed that our conclusions in this paper are not seriously 
affected even if we change the value of this parameter into 
a = 2.1 and 2.3. Then, the gamma-ray luminosity L^, which 
is defined as i/L^ in erg s^^ at 100 MeV (at the restframe), is 
related to the observed photon flux, at E-^ > =100 
MeV (in photons cm^^ s^' and Ey^i^ in the observer's frame), 
as 



where di is the standard luminosity distance. Then, the ob- 
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served data that will be used in the statistical analysis is a set 
of {zi,L^,i, Hi), where i7, denotes the observed blazar location 
in the sky, and the subscript / denotes each blazar, running 
over 1 < ! < A'obs = 46. 

2.2. Maximum Likelihood Method 

We constrain the model parameters of the blazar GLF mod- 
els by using the maximum likelihood method. A specified 
GLF model predicts the distribution function of the three 
quantities, d^N /{dz dL^ dil), and the likelihood function ^ 
is given as (see, e.g., Loredo & Lamb 1989): 



^-exp( yVexpj n ^^^^^^^ , 
where A^exp is the expected number of blazar detections: 



exp 



dz / dL^ / dil 



d^N 



dz dLj dfl 



(2) 



(3) 



Consider a transformation about the normalization, 
d'^N /{dz dLj dil) fd^N /{dz dL^ dVL). By maximizing 
the likelihood function about /, we find / = A^obs/A^exp, i-C-, 
the maximum likelihood obtained when the expected number, 
fNexp, becomes equal to A^obs- Substituting / — A^obs/-^exp 
and ignoring constant factors that are not relevant for the 
likelihood maximization, we find the normalization-free form 
of the likelihood function: 



Nobs 

n 

!=1 



1 d^N{zi,L^,i,ni] 



exp 



dz dL^ dil 



(4) 



(5) 



The distribution function can be expressed as 
d^N dV 

xe[F^{L^,z)~F^,iirn{^)] , 

where is the GLF per unit comoving density and unit lumi- 
nosity, dV /dz is the comoving volume element per unit solid 
angle as defined in the standard cosmology, and O is the step 
function [Q{x) — 1 and for x > and < 0, respectively]. 
The detection efficiency £{L^,z) represents the identification 
probability as a blazar by finding a radio counterpart, which 
will be defined in the next subsection. 

To find the best-fit model parameters and their confidence 
regions, we use the standard likelihood ratio method, assum- 
ing that J^f exp(— x^/2), where obeys the chi-square dis- 
tribution. The best-fit parameters are simply obtained as those 
giving the minimum chi-square, Xmin' confidence 
region is determined by the contour of a constant Ax^ = 
X^ — Xmin- ^'^ ^^^^ paper we will perform two-parameter fit to 
the data, and hence Ax^ obeys to the chi-square distribution 
with two degrees of freedom, i.e., Ax^ = 2.30, 6.16, and 9.21 
for 68%, 95%, and 99% C.L., respectively (see, e.g.. Press et 
al. 1992). 

2.3. Identification Probability in Radio Band 

Here we formulate the probability that a gamma-ray blazar 
having a gamma-ray flux above the EGRET detection limit 
is also detected in the radio band, so that it is identified as 
an EGRET blazar In the left panel of Figure |2l we show the 
observed relation between and radio luminosity Lr {vL^ 
in erg s^' at restframe 2.7 GHz) of the EGRET blazars. The 



best-fit relation is also shown in the figure {solid line). Here 
we assumed the photon spectral index of a, = 1-0 for the 
K-correction, as a typical index of blazars in the radio band 
(Miicke et al. 1997). It should be noted that, though the 
gamma-ray and radio luminosities are apparently well corre- 
lated with each other, this is mostly an artifact, coming from 
the fact that blazars with different distances are detected with 
a similar flux around the detection limit. This can easily be 
understood if we see the correlation plot between observed 
gamma-ray and radio fluxes, as shown in the right panel of 
the same figure. 

The correlation between gamma-ray and radio emissions of 
blazars has been investigated and the evidence for this cor- 
relation has been presented in many papers (e.g., Padovani 
et al. 1993; Stecker, Salamon, & Malkan 1993; Salamon & 
Stecker 1994; Dondi & GhiselHni 1995; Liihteenmaki et al. 
1997; Lahteenmaki, Valtaoja, & Tornikoski 2000; Tornikoski 
& Lahteenmaki 2000; Lahteenmaki & Valtaoja 2001). On 
the other hand, Miicke et al. (1997) claimed that correlation 
between the gamma-ray and radio luminosities cannot be es- 
tablished firmly from the existing data. This is probably be- 
cause the correlation is hidden by intrinsic dispersion and the 
rather narrow dynamic range of observed radio and gamma- 
ray fluxes. Based on this result, CM98 assumed no correlation 
between and Ly. 

However, it should be noted that the dynamic range of lu- 
minosity of EGRET blazars is extending over more than five 
order of magnitudes (Fig. |2}. If you assume no correlation 
between gamma-ray and radio luminosities, it means that you 
cannot tell which blazar is brighter in the radio band, even if 
you akeady know that one blazar is brighter than the other 
by a factor of 10^. Such an assumption is highly unlikely, 
since it is generally believed that the overall spectra of blazars 
are made by two different emission processes from the same 
population of relativistic electrons; the gamma-ray emission 
is due to the inverse Compton process, while the radio emis- 
sion is due to the synchrotron process. Therefore, we must 
introduce some correlation between the gamma-ray and radio 
luminosities in the analysis. 

Hence we introduce a linear correlation with log-normal 
scatter as a simple and phenomenological model, to avoid 
theoretical uncertainties of more physically motivated mod- 
els. Then, L^/Ly obeys to the log-normal distribution with 
(p) ^ 3.23, where p = \Q'gyQ{L^/Lr). Figure |3l shows the dis- 
tribution of p, and the best-fit dispersion is Op = 0.49. Then, 
the probability that a blazar having gamma-ray luminosity 
at redshift z will be identified in radio band can be calcu- 
lated as the probability that the corresponding radio flux at 
2.7 GHz (at observer's frame) is greater than the radio detec- 
tion limit, Frjim- We take /vum — 0.7 Jy, since most of the 
EGRET blazars have radio fluxes larger than 0.7 Jy. 

It should be noted that this correlation may be different 
from the true correlation between L-^ and L^, since the ob- 
served correlation has been affected by selection effects. The 
ratio of the flux limits in radio and gamma-ray bands (0.7 Jy 
at 2.7 GHz and - 1.0 x 10"^ photons cm"^ ^-i ^ihowe 100 
MeV, respectively) is in fact very close to the mean of p. In 
order to check how much our analysis could be affected by 
this effect, we calculated the prediction of the p distribution 
that will actually be observed for the EGRET blazars, from 
the best-fit models of the blazar GLF that will be obtained 
later in this paper. We confirmed that the predicted distribu- 
tion is consistent with the observed one, and this consistency 



4 



Narumoto and Totani 



check demonstrates that our analysis is not seriously biased 
by the selection effect. 

Though we assumed a simple linear relation, more physi- 
cally motivated models such as the synchrotron self-Compton 
(SSC) model and the external radiation Compton (ERC) 
model may predict deviation from the exact linear correlation. 
However, the above result indicates that the linear correlation 
is statistically consistent with the data, and we cannot derive 
more detailed conclusions from the current sample. The fu- 
ture GLAST mission may provide better statistics for this is- 
sue. We will discuss about the SSC and ERC models in the 
context of the beaming factor difference in the gamma-ray and 
radio bands in ^5.11 

2.4. Background Photon Flux from Unresolved Blazars 

We can calculate the integrated background photon flux (> 
100 MeV) from blazars below the EGRET detection limit as 



-^diffuse 



dz II. 



"y.lim 



(2) 



dL^ F^{L^,z) p^{L^,z) , (6) 



where L^.\im{z) is the gamma-ray luminosity corresponding 
to the EGRET threshold, and L^.min is the minimum gamma- 
ray luminosity of the blazar GLF. This quantity will be com- 
pared with the observed EGRB, to estimate the contribution 
from unresolved blazars. Since the minimum gamma-ray lu- 
minosity is quite uncertain and has considerable effect on the 
blazar contribution to the EGRB, we consider four values of 



10"*^ lO'*^ lO'*', and lO"*" erg s"'. For reference, 
10« 



^7,min 

io« erg s is smaller than the minimum gamma- 
ray luminosity of the EGRET blazars by a factor of 5. We 
assume Zrnax — 5, but the predicted EGRB flux hardly depends 
on this parameter, since the number density of AGNs with 
a given luminosity decreases with redshift beyond z ^ 2 by 
the estimated evolution of luminosity functions in the radio or 
X-ray bands, based on which our GLF models will be con- 
structed. 

3. MODELS OF THE BLAZAR GAMMA-RAY LUMINOSITY 
FUNCTION AND RESULTS OF THE ANALYSIS 

In this paper we try two models of the blazar GLF, 
Pj(Lj,z). The descriptions of these two models and fits to 
the observed data will be presented below. 

3.1. The Pure Luminosity Evolution (PLE) Model 

3.1.1. Model Description 

For the PLE model, we follow the same procedure pro- 
posed by SS96 for constructing the blazar GLF model. They 
made the basic assumptions that blazars seen by gamma-rays 
above 100 MeV are also seen in radio as FSRQs, and that the 
gamma-ray and radio luminosities of these objects are linearly 
related as 

= lO^'P^L,- , (7) 

where the definitions and units are the same as defined in the 
previous section. The blazar GLF is then derived from the 
FSRQ RLE: 



p^{L 



(8) 



where is a normalization factor, and pr{Lr,z) is the FSRQ 
RLE We use the FSRQ RLF derived by Dunlop & Peacock 
(1990, hereafter DP90): 



P,.(L,-,z) = ^P,.(^,0 



where /9,-(L, ,0) is the present-day FSRQ RLF, which is char- 
acterized by the faint-end slope index 71, the bright-end slope 
index 72, and the break luminosity L*, given as 



PriLr,0) = 



Ar 



r 




71 




72 ^ 






+ 







(In 10) L, 

and f{z) is the luminosity evolution function given as 

f{z) - 10"^+''^' . 



(10) 



(11) 



(9) 



Here, A, = 7.08 x 10"^ Mpc-^ log^L; 42.79, 71 = 0.83, 
72 = 1.96, a = 1.18, and b = -0.28. Since this FSRQ RLF 
was derived for the Einstein-de Sitter (EdS) universe with 
(17o,fiA) = (1.0,0.0) and//o = 50kms-i Mpc"^ we multi- 
ply pr{Lr,z) by a correction factor {dVEds/dVA) and f{z) by 
a correction factor {di.K/dL.Ed^Y in order to approximately 
convert to a FSRQ RLF for the ACDM universe, where dV is 
the comoving volume element of the universe, and di is the 
luminosity distance. Therefore, this model is no longer "PLE" 
model in a strict sense, but we take this correction into account 
to apply the RLF in agreement with the observed data. 

3.1.2. Constraints from the Redshift and Luminosity 
Distribution of the EGRET blazars 

In this model, we take {p) and 71 as the two free param- 
eters since they are poorly constrained from observations, 
and fix the other parameters to the best-fit values shown in 
^3.1.11 For consistency, we use the same {p) with the disper- 
sion dp = 0.49 also for judge ment of the radio identification, 
which has been described in ^2.31 The normalization factor 
77, which is physically related to the possible beaming effect, 
is determined by the requirement that the calculated number 
of identifiable blazars above the EGRET threshold is equiva- 
lent to the observed number of the EGRET blazars. In Fig- 
urelUwe show the 68%, 95%, and 99% C.L. contours for the 
PLE model parameters {solid lines). The best-fit parameters, 
i{p)^li) = (3.28,0.69), are also marked {cross). The value 
{p) = 3.28 is quite similar to the value directly obtained from 
the EGRET blazars {p — 3.23). The faint end slope 71 — 0.69 
is somewhat smaller (i.e., flatter faint-end slope) than that of 
the FSRQ RLF derived by DP90 (71 = 0.83), but the value of 
DP90 is well within the 68% C.L. region. 

Figures |5] and |6l show the redshift and luminosity distribu- 
tions for the best-fit parameters, respectively {dashed lines). It 
is clear that the PLE model with parameters adopted by SS96 
{{p) = 2.54 and 71 = 0.83 from DP90) can reproduce neither 
the redshift nor luminosity distributions. Our best-fit model 
reproduces these distributions better than the SS96 model, but 
still the fit is not very good, especially for the redshift distribu- 
tion. We performed the Kolmogorov-Smirnov (KS) test, and 
find that the chance probability of getting the observed devi- 
ation of the redshift distribution from the best-fit PLE model 
is only 3.1%, while it is 27.0% for the luminosity distribu- 
tion. These results indicate that the PLE framework may not 
be satisfactory to describe the EGRET blazar data. 

3.1.3. Blazar Contribution to the Extragalactic Diffuse 
Gamma-Ray Background 

In Figure we present the contours of 25%, 50%, 75%, 
and 100% contribution of blazars under the EGRET sensitiv- 
ity limit to the EGRB for the PLE model {dashed lines). We 
find that unresolved blazars can explain only 50-55% of the 
EGRB for the best-fit parameters. On the other hand, since 
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the contour of 100% blazar contribution pass through inside 
the 68% C.L. region for all the cases, it is unable to exclude 
the possibility that almost all of the EGRB is explained by 
blazars. However, the poor fit of the PLE model to the ob- 
served redshift distribution indicates that it is not appropriate 
to derive any conclusion about the EGRB based on this model 
framework. It is apparent from this figure that the blazar con- 
tribution to the EGRB is strongly dependent on the faint end 
slope 71. On the other hand, since the EGRB contribution 
becomes 100% in a region where 71 < 1.0, the minimum 
gamma-ray luminosity of blazars, L^.min. hardly affects the 
contribution to the EGRB. 

3.2. The Luminosity Dependent Density Evolution ( LDDE) 

Model 

3.2.1. Model Description 

In this section, we construct the blazar GLF model based 
on the AGN SXLF by assuming a linear relation between the 
blazar gamma-ray luminosity (dominated by the jet) and the 
AGN soft X-ray luminosity (dominated by the disk emission) 
expressed as 

= lO^Lx , (12) 
where the unit of (in lyL^, at the restframe 100 MeV) and 
Lx (in the restframe 0.5-2 keV X-ray band) is erg s^' . In the 
soft X-ray bands, the typical AGN spectra have a photon in- 
dex of ~ 2, i.e., constant in lyF^, and hence Lx in the observed 
0.5-2 keV band is used as that in the restframe 0.5-2 keV 
band (e.g., Hasinger, Miyaji & Schmidt 2005). The assump- 
tion of the linear relation between Lx and L~^ is motivated 
from an expectation that the jet activity should be somehow 
correlated with the accretion power which can be measured 
by X-ray luminosity from accretion disks. It should be noted 
that this relation is not necessarily be seen in the observed 
spectral energy distributions of blazars, since X-ray emission 
from blazars is dominated by the beamed emission from the 
jet, rather than the disk emission. 

The blazar GLF is then obtained from the AGN SXLF: 



(13) 



where k is a normalization factor, and px{Lx,z) is the AGN 
SXLF. In this model, we adopt the same form as derived by 
Hasinger, Miyaji, & Schmidt (2005) for the AGN SXLF and 
the details are as follows: 

px{Lx,z)^PxiLx,0)fiLx,z) , (14) 

where px {Lx , 0) is the present-day AGN SXLF, which is char- 
acterized by the faint-end slope index 71, the bright-end slope 
index 72, and the break luminosity L^, given as 

and f{Lx,z) is the density evolution function given as 



(1 



nLx,z) = ( f^LxMLx) 



1+z 



l+ZciLx) 



P2 



[z<ZciLx)] , 
[z>ZciLx)] , 
(16) 



where Zc is the redshift of evolutionary peak given as 

Z* {Lx>La), 

z* 



Zc{Lx) 



and p\, P2 given as 

p\ =Pi 



-A(logioix-44), 
-/32(logioix-44). 



(18) 
(19) 



Here, the parameters obtained by the fit to X-ray data are 
(Hasinger, Miyaji, & Schmidt 2005): Ax = 6.69 x 10"^ 
Mpc-\ logioL^. = 43.94 ±0.11, 71 = 0.87 ±0.10, 72 = 
2.57±0.16, z* = 1.96±0.15, logioi„ = 44.67, a = 0.21 ± 
0.04, pI =4.7±0.3, pI = -1.5±0.7, A =0.7±0.3, and 
/32 = 0.6 ±0.8. 

3.2.2. Constraints from the Redshift and Luminosity 
Distribution of the EGRET blazars 

In this model, we take q and 71 as the two free parame- 
ters and fix the rest to the best-fit parameters described in the 
previous section. For the radio identification judgement, we 
use the value of {p) = 3.23 as obtained from the L^ — L,- re- 
lation of the EGRET blazars. The normalization factor k is 
determined by fitting the expected total number of blazars to 
the actually observed number of the EGRET blazars. In Fig- 
ure El we show the 68%, 95%, and 99% C.L. contours for 
the LDDE model {solid lines), with the best-fit parameters 
(^,71) = (3.80, 1.19) marked by cross. The best-fit value of 
7i is slightly larger than the value inferred from the SXLF 
(71 = 0.87 ± 0.10), but the SXLF value is within the arrowed 
region of ^ 95% confidence level. Figures [S] and |6l show the 
redshift and luminosity distributions for the best-fit parame- 
ters, respectively (solid lines). It is noteworthy that the LDDE 
model can reproduce the redshift and luminosity distributions 
of the EGRET blazars better than the PLE model. Quantita- 
tively, the chance probability of getting the observed deviation 
estimated from the KS test is 67.8% and 99.3% for the redshift 
and luminosity distributions, while these are 3.1% and 27.0% 
for the best-fit PLE model, respectively. These results indi- 
cate that the blazar evolution can better be described by the 
LDDE rather than the PLE. 

3.2.3. Blazar Contribution to the Extragalactic Diffuse 
Gamma-Ray Background 

In Figure Owe present the contours of 25%, 50%, 75%, 
and 100% blazai- contribution to the EGRB for the LDDE 
model {dashed lines). We find that unresolved blazars can ex- 
plain only 25-50% of the EGRB for the best-fit parameters. 
For the case of L-y ^in — lO'*^ erg s^', the contour of 100% 
blazar contribution is outside the 99% C.L. region. Still, if we 
take the case of L^,mi„ = lO'*" erg s"', the LDDE GLF with 
71 ^ 1.26 can marginally explain 100% of the EGRB with 
the parameters within the 68% C.L. region. However, such 
a steep faint-end slope index is not favored from the SXLF 
(71 = 0.87 ± 0.10). In the LDDE model the blazar contri- 
bution to the EGRB is strongly dependent on the minimum 
gamma-ray luminosity of blazars L^_min as well as the faint 
end slope 71 , since the best-fit faint end slope is 71 > 1 .0. 

4. PREDICTIONS FOR THE GLAST MISSION 

4.1. Expected Number of GLAST Blazars and Their 
Contribution to the EGRB 

The number of blazars with flux stronger than can be 
calculated as 



dV f°° 

^^~r I Pii^n^) ' (20) 

dz Jl^{z.,f^) 
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where L^{z,F^) is the gamma-ray luminosity of a blazar at 
redshift z whose photon flux above 100 MeV is F-^ (see eq. 
(!]). In the left panel of Figure^! we show the calculated 
logA^ — log/^y relation of blazars. This figure shows that the 
SS96, our best-fit PLE, and the LDDE models predict consid- 
erably different numbers of blazars detectable by the GLAST 
(~ 10000, 5350 and 3000, respectively), where we have set 
the GLAST sensitivity limit as — 2.0 x 10^^ photons 
cm^^ s^^ Here we used L-y.min = 10^" erg s^', but the de- 
pendence of the predicted counts on this parameter is small. It 
is remarkable that the LDDE model predicts more than three 
times fewer blazars than the SS96 model. This is because 
the LDDE model predicts smaller evolution for less luminous 
blazars, and hence paucity of high-z and faint blazars, which 
have the dominant contribution to the blazar counts at faint 
flux in the SS96 or our best-fit PLE models. This means that 
we can constrain different blazar GLF models and their cos- 
mological evolution from the number counts of blazars de- 
tected by the GLAST, even without knowing their redshifts. 
We also calculate the predicted counts with the LDDE model 
parameters of (^,71) = (3.80, 1.26), which are within the ar- 
rowed region of the 68% C.L. and able to explain 100% of 
the EGRB. In this case the prediction for the GLAST is in- 
creased to ~ 4700, but still smaller than those of the SS96 or 
the best-fit PLE model. 

How much fraction of the EGRB can be resolved by the 
GLAST mission? To answer to this question, in the right 
panel of Figure^] we show the differential flux distribution 
of gamma-ray blazars multiplied by flux, showing the contri- 
bution to the EGRB per unit logarithmic flux interval. For the 
PLE model, we predict that we will see the peak of the con- 
tribution to the EGRB above the detection limit of the GLAST 
mission, and hence we will resolve a considerable fraction 
of the EGRB into blazars, if unresolved blazars are the ma- 
jor source of the EGRB. The predicted resolvable fraction 
of the EGRB flux by blazars detectable by the GLAST (but 
under the EGRET detection Umit) is 33% and 42% for the 
best-fit PLE model and that with ((/?}, 71) = (3.28,0.85), re- 
spectively. The latter model can explain 100% of the EGRB 
by unresolved blazars. On the ohter hand, the LDDE model 
curves have two peaks of the contribution to the EGRB as 
a function of F^, because of the complicated nature of the 
evolution. We predict that the contribution to the EGRB will 
decrease with decreasing flux, just below the EGRET sen- 
sitivity limit. The resolvable fraction of the EGRB by the 
GLAST is 20% and 26% for the best-fit LDDE model and that 
with {q, 71 ) = (3.80, 1 .26), where the latter model can explain 
100% of the EGRB. As shown in Figure[^ the dominant con- 
tribution to the EGRB comes from blazars under the GLAST 
detection limit, even if blazars are the dominant source of the 
EGRB. 

4.2. Redshift and Luminosity Distribution 

In Figures [TTI and [T2I we show the redshift and luminosity 
distributions of blazars detectable by the GLAST, respectively. 
It should be noted that only the shapes of distribution should 
be compared, since the total number has been normalized to 
the same. It can be seen that the peak of the redshift distri- 
bution in the LDDE model occurs at a lower redshift than the 
best-fit PLE or SS96 models. Furthermore, both redshift and 
luminosity distributions of the LDDE model are wider than 
those of the other two models. Though the redshift must be 
determined for the future GLAST blazars, this will provide 
another important information to discriminate different GLF 



models. 



5. DISCUSSION 



5.1. Normalization, Beaming, and Duty Cycle 

In Figures |3 and |8l we also present the contours of 77 and 
K, {dashed lines). These parameters are possibly related to 
a beaming effect and/or duty cycle, and we find r/ ^ 10^"^ 
and K ~ 10^^'^ for the best-fit PLE and LDDE models, re- 
spectively. The inferred value of 7/ is not far from the unity, 
as expected because the GLF was constructed from the lu- 
minosity function of FSRQs, which are generally believed to 
be the same population with blazars. However, the best-fit 
value -q ~ lO^*^'^ is slightly smaller than the unity. In the 
SSC model, this value may be explained if not all FSRQs 
are sufficiently gamma-ray-loud, or there are more than one 
components of nonthermal electrons or emission regions hav- 
ing different beaming patterns (e.g., Lindfors et al. 2005). 
On the other hand, the ERC model may explain this value 
of 7] by a single electron component, since the beaming pat- 
tern of the ERC emission is narrower than that of the syn- 
chrotron emission or the SSC emission (Dermer 1995). The 
observed flux has an angular dependence Su °^ for the 

synchrotron or SSC processes, and Si, °^ f^2+2a ^^j. ^j^g ^^^q 
process, where ^ is the Doppler factor and a is the photon 
spectral index. Here we define 9e as the viewing angle of the 
observer measured from the jet axis, at which the observed 
flux is smaller than that for the direction of the jet axis by a 
factor of e. Then, using the Lorentz factor F = 10 and the 
typical index (a^ = 2.2, = 1.0), we find that 9^ ^ 3.6° for 
the synchrotron emission, while 9^ ^2.1° for the ERC emis- 
sion. Therefore 77 ^ (2. 1 /3.6)^ ~ 0.34 is expected in this case, 
which is moderately close to 77 ~ 10^"^. 

On the other hand, the GLF in the LDDE model is con- 
structed from that of AGNs selected in the soft X-ray band, 
most of which are expected to be type 1 AGNs. Since the jet 
activity is not always observed in AGNs, it is rather unlikely 
that all X-ray AGNs have blazar activity. Then, the interpreta- 
tion of the small value of k is a combination of the duty cycle 
^ (here defined as the fraction of AGNs having the blazar ac- 
tivity) and beaming of radiation, i.e.. 



/typel 



= 5 X 10" 



An 

47r 

-6 I /typel 

0.2 



/AO/(47r) 



10 



-3 



10 



-3 



(21) 
(22) 



Here, we take the ratio of SXLF to HXLF normalization as the 
fraction of type 1 AGNs in all AGNs, /typei (Hasinger, Miyaji, 
& Schmidt 2005). The beaming expected from the estimated 
Lorentz factor of blazars (F ~10-20, e.g., Maraschi & Tevec- 
chio 2003) is Arj/(47r) - l/(4r2) - 10"^ The smafl duty 
cycle can partially be ascribed to the fraction of radio-loud 
AGNs to all AGNs, /radio ~ 0.15 (Urry & Padovani 1995; 
Krolik 1999). Comparing this /radio to the inferred ^, it is in- 
dicated that about 1% of radio-loud galaxies have active jets 
now. It is not unreasonable, since the jet activity may be spo- 
radic, and radio emission by relic electrons can be kept for a 
while after the jet activity ceased. 

5.2. Unidentified EGRET Sources 

Over half of the gamma-ray sources (170 of 271) detected 
by the EGRET have not been identified as known astronom- 
ical objects. The distribution of these unidentified sources 
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can be accounted for as the sum of the Galactic and an- 
other isotropic (Hkely extragalactic) component (Mukherjee 
et al. 1995; Ozel & Thompson 1996). Since most of the 
firmly identified extragalactic sources are blazars, it is widely 
believed that many of the unidentified extragalactic sources 
which show high variability are blazars. Actually, many re- 
searchers have been investigating the unidentified EGRET 
sources in various wavelengths to find their counterparts, and 
some of them, such as 3EG J2016H-3657 (Mukherjee et al. 
2000; Halpem et al. 2001), BEG J2006-2321 (Wallace et al. 
2002), and BEG J2027-I-3429 (Sguera et al. 2004), have al- 
ready been identified as blazars. In addition, the number of 
blazar candidates in the unidentified sources is now increasing 
(e.g., Sowards-Emmerd, Romani, & Michelson 2003; Bloom 
et al. 2004; Wallace, Bloom, & Lewis 2005). 

These (potential) blazars were initially classified as uniden- 
tified sources likely because of the lack of strong radio emis- 
sion. Since our model incorporates the dispersion in radio 
and gamma-ray luminosities and selection by radio flux, we 
can address this issue by estimating how many "unidentified 
blazars" will be expected due to the lack of detectable ra- 
dio flux. We found that this number is 10 and 8 sources at 
high Galactic latitude of > 45° for the best-fit PLE and 
LDDE models, respectively^. In the third EGRET catalog, 
there are six low confidence potential blazars and 19 uniden- 
tified sources at \b\ > 45°. The variability of these uniden- 
tified sources have been investigated in some papers (e.g., 
Gehrels et al. 2000; Torres, Pessah, & Romero 2001; Nolan et 
al. 2003), and 5-9 of them are non-variable sources, though 
definition of "non-variable" is different among these papers. 
Therefore, a substantial fraction of potential blazars and vari- 
able unidentified sources at \b\ > 45° can be explained by 
blazars with radio flux under the detection hmit. It has been 
suggested that the apparently steady unidentified sources may 
be accounted for by forming gamma-ray clusters (Totani & 
Kitayama 2000). 

6. CONCLUSIONS 

In this paper, we presented a comprehensive study for the 
gamma-ray luminosity function of blazars. We introduced a 
log-normal distribution for the ratio of radio to gamma-ray 
luminosity, and radio detection was required for the model 
blazars to become identified sources in the EGRET catalog. 
The number of potential blazars and variable unidentified 
sources at high Galactic latitude in the EGRET catalog is sim- 
ilar to the number of "unidentified blazars" predicted in our 
model due to the lack of detectable radio flux, indicating that 
our treatment about gamma-ray blazar identification by radio 
detection is reasonable. We newly tried the luminosity depen- 
dent density evolution (LDDE) model based on recent studies 
of the X-ray luminosity function of AGNs, in addition to the 
pure luminosity evolution (PLE) model used in earlier studies. 

By performing the maximum likelihood analysis for the 
redshift and luminosity distributions of the EGRET blazars, 
we found that the LDDE model with the evolutionary param- 
eters inferred from the soft X-ray luminosity function (SXLF) 
of AGNs can explain the redshift and luminosity distributions 
of the EGRET blazars better than the PLE model with the evo- 
lutionary parameters inferred from the radio luminosity func- 
tion (RLE) of flat-spectrum radio-loud quasars (FSRQs). This 
indicates that blazars are evolving similarly to type 1 AGNs 

^ The Galactic component of unidentified sources extends only to |6j ~ 45° 
(Gehrels et al. 2000). 



found in the soft X-ray bands, and hence the jet activity is uni- 
versally correlated with the accretion history of AGNs. We 
also found that the normaUzation between blazars and type 1 
AGNs is roughly consistent with the unified picture of AGNs, 
when the beaming and jet duty cycle are taken into account. 

As an implication for the GLAST mission, we found that the 
LDDE model predicts considerably fewer (by a factor of more 
than 3) blazars down to the GLAST sensitivity limit, com- 
pared with a previous estimate based on the PLE luminosity 
function. This can be easily tested by the mission, giving us 
important information for the evolutionary nature of gamma- 
ray blazars. Redshift and luminosity distributions will further 
constrain the different models, though redshift measurements 
are required. 

Then we examined the contribution of unresolved blazars 
to the extragalactic diffuse gamma-ray background (EGRB), 
which has been controversial topic in earlier works (100% in 
SS96, 25% in CM98). We found that only 25-50% of the 
EGRB can be explained with the best-fit LDDE model, which 
is similar to the result of CM98 but by considerably different 
analysis. On the other hand, according to our statistical analy- 
sis and parameter survey, it is possible to account for 100% of 
the EGRB with a steeper faint-end slope of 71 ~ 1 .26 that is 
marginally consistent with the arrowed region from the Uke- 
Uhood analysis of the EGRET blazar distributions. However, 
such a value is inconsistent with that inferred from the SXLF 
of AGNs. Therefore we conclude that unresolved blazars can- 
not account for 100% of the EGRB, if the jet activity of AGNs 
is universally correlated to the accretion luminosity and hence 
the AGN SXLF is a good description of the blazar luminosity 
function and its evolution. It should be noted that the uncer- 
tainty about the extrapolation to high redshift does not change 
this conclusion, since almost all of the cosmic X-ray back- 
ground (CXB) flux can be explained by the LDDE luminosity 
function (Ueda et al. 2003). It indicates that, if the rest of the 
EGRB is explained by an AGN population, it must be a dif- 
ferent population from EGRET blazars, having different evo- 
lution from X-ray AGNs, and not significantly contributing to 
flieCXB. 

Based on the best-fit LDDE model, we predict that the con- 
tribution to the EGRB by blazars will start to decrease with 
decreasing flux just below the EGRET sensitivity limit. In 
the case of the LDDE model with parameters that can explain 
100% of the EGRB, there are two peaks of the contribution 
to the EGRB as a function of flux, and the major contribution 
comes from blazars under the GLAST detection Umit. There- 
fore it is unlikely that almost all the EGRB flux is resolved 
into discrete blazars even by the sensitivity of GLAST. 

This work was supported by the Grant-in-Aid for the 
21st Century COE "Center for Diversity and Universality in 
Physics" from the Ministry of Education, Culture, Sports, Sci- 
ence and Technology (MEXT) of Japan. 
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Fig. 1 . — Left panel: Sensitivity limit to point sources as a function of the Galactic latitude in the third EGRET catalog. The fitted relation is also shown by the 
sohd line. Right panel: mean flux of the EGRET blazars against its Galactic latitude. The solid fine is the same as the left panel. 
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Fig. 2. — Left (Right) panel: Observed gamma-ray and radio luminosities (fluxes) of the EGRET blazars are shown by the crosses. Here, the luminosity (flux) 
is defined as vL^ (vF^,) at restframe (observed) 100 MeV and 2.7 GHz for the gamma-ray and radio bands, respectively. The solid fine is the best-fit linear 
relation. The K-corrections are done for luminosities assuming typical spectral indices (see text). 
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Fig. 3. — Histogram of radio to gamma-ray luminosity ratio, p = log,iQ{L^/Lr), of the EGRET blazars. The luminosities are vL^ in the restframe 100 MeV 
and 2.7 GHz bands, respectively. The soUd curve is a Gaussian fit to the histogram. 
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Fig. 4. — The solid contours show the 68%, 95% and 99% C.L. regions for the PLE model parameters [the faint-end slope index 71 and the mean gamma-ray 
to radio luminosity ratio, (p) = (log[g(L-,/Lr))]. The best-fit values, ((p),7i) = (3.28,0.69), are shown by the cross. The dashed contours correspond to 
T] = 10""'^^, 10"'^-**, lO^'-", and 10"'-'^, respectively, where r] is the ratio of the normalizations of the gamma-ray to radio luminosity functions. 
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Fig. 5. — Redshift distribution of the EGRET blazars. The histogram is the EGRET data. The solid and dashed curves are the best-fit models for the LDDE and 
PLE models, respectively, from the likeUhood analysis. The dotted curve is obtained from the blazar GLF model of SS96. The error bars are \a Poisson error. 
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Fig. 6. — Luminosity distribution of tlie EGRET blazars. The line markings are tlie same as Figure|5] The luminosity is uL^, at 100 MeV. The en'or bars are 
la Poisson error 
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Fig. 8. — The solid contours show the 68%, 95% and 99% C.L. hkehhood contours for the LDDE model parameters [the faint-end slope index 71 and the 
gamma-ray to X-ray luminosity ratio, q = log[o(L-,/Lx)]- The best-fit values, (g,7i ) = (3.80, 1.19), are shown by the cross. The dashed contours correspond 
to K = 10^'*", 10^**^^, lO^'*'**, 10^^", 10"'^^^, 10"'^*, and 10^*", respectively, where k is the normalization ratio of the gamma-ray to soft X-ray luminosity 
functions. 
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Fig. 9. — The solid contours and crosses are the same as Figure|^ showing the fit by the LDDE model. The dashed contours show 25%, 50%, 75%, and 100% 
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Fig. 10. — Left panel: logA? — logF distribution (cumulative flux distribution) of blazars. The solid and dashed curves show the prediction by the best-fit LDDE 
and PLE models, respectively. The dotted curve is derived from the blazar GLF model of SS96. The observed distribution of the EGRET blazars is shown by the 
thin solid line. The detection limits of the EGRET and GLAST are also shown in the figure. Right panel: The same as the left panel, but showing differential flux 
distribution multiplied by F-y , to show the contribution to the EGRB per logarithmic flux interval. The thick solid and dashed curves are the same as those in the 
left panel, but the thin solid and dashed curves show the LDDE and PLE models with parameters that can explain all the EGRB flux by unresolved blazars. (See 
the labels in the panel for the values of the parameters.) 
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Fig . 11. — Redshift distribution of blazars expected for the GLAST observation. The total number is normalized to the same. The solid and dashed curves show 
the prediction by the best-fit LDDE and PLE models, respectively. The dotted curve is predicted from the blazar GLF model of SS96. The GLAST sensitivity 
limit is set as Fy„^ = 2.0 x 10"' photons cm"^ s"'. 
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Fig. 12. — Luminosity distribution of blazars expected for the GLAST observation. The total number is normalized to the same. The line markings are the 
same as Figure lTTI The GLAST sensitivity limit is set as Fii„, = 2.0 x 10"' photons cm"^ s"' . 



